clear 

disp('Baseline')

%% Parameter choices

gamma = 3;                      % risk-aversion coefficient
    u = @(x) x.^(1-gamma)/(1-gamma);    % utility function: assume CRRA utility based on terminal wealth.
T = 40;                         % years between first contribution and retirement.
I_0 = 1;                        % PH's initial contribution to QRP.
g_I = 0.03;                     % annual growth rate of investor's contribution.

fm_setting = "baseline";        % financial market setting ("baseline", "high interest rate", "pessimistic").

% RILA characteristics:
tau = 1;                        % length of RILA crediting term (in years).
RILA_fee_rate = 0;              % RILA fee rate.
c_bar = 0;                      % expected cost (p.a.) of RILA contract. That is: PH gets 100*(1-c_bar)\% of her investment back (in risk-neutral PV terms).

% Mutual fund:
MF_fee_rate = 0;                % Mutual fund (LSF / TDF) fee rate.


%% Obtain Gauss-Hermite Quadrature nodes & weights for fund returns
switch fm_setting
    case "baseline"
        r = 0.03;
        lambda = 68.150;
        mu_d_RN = 0.1629 + r;
        MRP = 0.055;
        sigma_d = 0.10343;
        mu_j = -0.0025643;
        sigma_j = 0.018482;
    case "high interest rate"
        r = 0.05;
        lambda = 68.150;
        mu_d_RN = 0.1629 + r;
        MRP = 0.055;
        sigma_d = 0.10343;
        mu_j = -0.0025643;
        sigma_j = 0.018482;
    case "pessimistic"
        r = 0.02;           % constant annual interest rate, compounded continuously.
        lambda = 0.29;      % likelihood of jump.
        mu_d_RN = 0.1148;
        MRP = 0.035;
        sigma_d = 0.10;     % volatility of diffusion process.
        mu_j = -0.4;        % mean of each jump.
        sigma_j = 0.09;     % volatility of magnitude of each jump.
    otherwise
        disp('Financial market setting unspecified.')
end

[R_S_GHQ,R_S_GHQ_RN,weight] = CreateNetReturns_GHQ(lambda, mu_d_RN, MRP, sigma_d, mu_j, sigma_j);

A_grid = I_0*[0:0.2:5,6:1:100,100*exp(0.02*(1:225))]';      % Dynamic programming grid of QRP account value


%% Find certainty equivalent of each investment option
tic

CE = zeros(11,1);

% [0] Risk-free investment.
Alpha_choices = 0;
[CE(1),~] = Get_CE_LSF_optim_fun(Alpha_choices,MF_fee_rate,u,T,g_I,R_S_GHQ,weight,r,A_grid,I_0);

% [1] Life-Style Fund with time-invariant (optim.) equity exposure 'Alpha'.
Alpha_choices = 0:0.01:1;
[CE(2),Alpha_star_LSF] = Get_CE_LSF_optim_fun(Alpha_choices,MF_fee_rate,u,T,g_I,R_S_GHQ,weight,r,A_grid,I_0);

% [2] RILA with time-invariant (optim.) Buffer 'B'.
Floor_control = 0; Buffer_control = 1;
F_choices = 1; B_choices = 0:0.01:0.50; Pu_choices = 1; Pd_choices = 1; Pu_Pd_same = 1;
[CE(3),~,B_star_fixed,~,~] = ...
    Get_CE_RILA_optim_fun(Floor_control,F_choices,Buffer_control,B_choices,Pu_choices,Pd_choices,Pu_Pd_same,u,RILA_fee_rate,c_bar,T,g_I,R_S_GHQ,R_S_GHQ_RN,weight,r,tau,A_grid,I_0);

% [3] RILA with time-invariant (optim.) Floor 'F'.
Floor_control = 1; Buffer_control = 0;
F_choices = 0:0.01:1; B_choices = 0; Pu_choices = 1; Pd_choices = 1; Pu_Pd_same = 1;
[CE(4),~,F_star_fixed,~,~] = ...
    Get_CE_RILA_optim_fun(Floor_control,F_choices,Buffer_control,B_choices,Pu_choices,Pd_choices,Pu_Pd_same,u,RILA_fee_rate,c_bar,T,g_I,R_S_GHQ,R_S_GHQ_RN,weight,r,tau,A_grid,I_0);

% [4] Target-Date Fund with dynamically optimized equity exposure 'Alpha'.
Alpha_choices = 0:0.01:1;
[CE(5),Alpha_star_mat_TDF,Alpha_star_0_TDF] = Get_CE_TDF_optim_fun(Alpha_choices,MF_fee_rate,u,T,g_I,R_S_GHQ,weight,r,A_grid,I_0);

% [5] Target-Date RILA with Buffer feature and dynamically optimized rates 'B'. 
Floor_control = 0; Buffer_control = 1;
F_choices = 1; B_choices = 0:0.01:0.50; Pu_choices = 1; Pd_choices = 1; Pu_Pd_same = 1;
[CE(6),~,B_star_mat_TDRILA,~,~,~,B_star_0_TDRILA,~,~] = ...
    Get_CE_TDRILA_optim_fun_parfor(Floor_control,F_choices,Buffer_control,B_choices,Pu_choices,Pd_choices,Pu_Pd_same,u,RILA_fee_rate,c_bar,T,g_I,R_S_GHQ,R_S_GHQ_RN,weight,r,tau,A_grid,I_0);

% [6] Target-Date RILA with Floor feature and dynamically optimized rates 'F'.  
Floor_control = 1; Buffer_control = 0;
F_choices = 0:0.01:1; B_choices = 0; Pu_choices = 1; Pd_choices = 1; Pu_Pd_same = 1;
[CE(7),F_star_mat_TDRILA,~,~,~,F_star_0_TDRILA,~,~,~] = ...
    Get_CE_TDRILA_optim_fun_parfor(Floor_control,F_choices,Buffer_control,B_choices,Pu_choices,Pd_choices,Pu_Pd_same,u,RILA_fee_rate,c_bar,T,g_I,R_S_GHQ,R_S_GHQ_RN,weight,r,tau,A_grid,I_0);

% [7] Enhanced Target-Date RILA with Buffer feature and dynamically optimized rates 'B', 'Pu' = 'Pd'. 
Floor_control = 0; Buffer_control = 1;
F_choices = 1; B_choices = 0:0.01:0.50; Pu_choices = 0.5:0.05:1; Pd_choices = 1; Pu_Pd_same = 1;
[CE(8),~,~,~,~,~,~,~,~] = ...
    Get_CE_TDRILA_optim_fun_parfor(Floor_control,F_choices,Buffer_control,B_choices,Pu_choices,Pd_choices,Pu_Pd_same,u,RILA_fee_rate,c_bar,T,g_I,R_S_GHQ,R_S_GHQ_RN,weight,r,tau,A_grid,I_0);

% [8] Enhanced Target-Date RILA with Floor feature and dynamically optimized rates 'F', 'Pu' = 'Pd'. 
Floor_control = 1; Buffer_control = 0;
F_choices = 0:0.01:1; B_choices = 0; Pu_choices = 0.5:0.05:1; Pd_choices = 1; Pu_Pd_same = 1;
[CE(9),~,~,~,~,~,~,~,~] = ...
    Get_CE_TDRILA_optim_fun_parfor(Floor_control,F_choices,Buffer_control,B_choices,Pu_choices,Pd_choices,Pu_Pd_same,u,RILA_fee_rate,c_bar,T,g_I,R_S_GHQ,R_S_GHQ_RN,weight,r,tau,A_grid,I_0);
 
% [9] Enhanced Target-Date RILA with Buffer feature and dynamically optimized rates 'B', 'Pu' and 'Pd', with leverage. 
Floor_control = 0; Buffer_control = 1;
F_choices = 1; B_choices = 0:0.01:0.50; Pu_choices = 0.5:0.05:1.5; Pd_choices = 0.5:0.05:1.5; Pu_Pd_same = 0;
[CE(10),~,B_star_mat_ETDRILA,Pu_star_mat_ETDRILA_B,Pd_star_mat_ETDRILA_B,~,B_star_0_ETDRILA,Pu_star_0_ETDRILA_B,Pd_star_0_ETDRILA_B] = ...
    Get_CE_TDRILA_optim_fun_parfor(Floor_control,F_choices,Buffer_control,B_choices,Pu_choices,Pd_choices,Pu_Pd_same,u,RILA_fee_rate,c_bar,T,g_I,R_S_GHQ,R_S_GHQ_RN,weight,r,tau,A_grid,I_0);

% [10] Enhanced Target-Date RILA with Floor feature and dynamically optimized rates 'F', 'Pu' and 'Pd', with leverage. 
Floor_control = 1; Buffer_control = 0;
F_choices = 0:0.01:1; B_choices = 0; Pu_choices = 0.5:0.05:1.5; Pd_choices = 0.5:0.05:1.5; Pu_Pd_same = 0;
[CE(11),F_star_mat_ETDRILA,~,Pu_star_mat_ETDRILA_F,Pd_star_mat_ETDRILA_F,F_star_0_ETDRILA,~,Pu_star_0_ETDRILA_F,Pd_star_0_ETDRILA_F] = ...
    Get_CE_TDRILA_optim_fun_parfor(Floor_control,F_choices,Buffer_control,B_choices,Pu_choices,Pd_choices,Pu_Pd_same,u,RILA_fee_rate,c_bar,T,g_I,R_S_GHQ,R_S_GHQ_RN,weight,r,tau,A_grid,I_0);

CE

toc

% save Results_baseline_noFees.mat